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►  We  succeed  in  reducing  computation  cost  without  losing  simulation  accuracy. 
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Improving  the  accuracy  of  the  sub-grid-scale  (SGS)  model  designed  for  three-dimensional  (3D) 
numerical  simulation  of  solid  oxide  fuel  cell  (SOFC)  electrodes  is  accomplished  by  introducing  an 
adaptive  power  index,  which  is  determined  by  the  length  ratio  between  the  grid  size  and  the  charac¬ 
teristic  scale  of  the  original  porous  structure.  The  main  function  of  the  SGS  model  is  to  maintain  the 
quality  of  the  voxel-based  structural  information  in  a  resampled  structure  used  as  the  calculation  grid 
system  by  considering  a  structure  with  a  characteristic  scale  that  is  smaller  than  the  grid  size. 
Improvement  of  the  simulation  accuracy  and  reduction  of  the  computational  load  can  be  expected  with 
the  SGS  models.  In  this  study,  four  types  of  SGS  models  are  investigated  to  clarify  their  effect  and 
applicability.  They  are  featured  by  (1)  volume  conservation,  (2)  interfacial  connectivity,  (3)  power  law 
with  a  constant  power  index,  and  (4)  power  law  with  an  adaptive  power  index,  respectively.  The  latter 
two  models  are  based  on  a  power-law  expression  of  the  local  tortuosity  factor  and  are  newly  proposed  in 
this  study.  Although  the  first  three  models  have  a  positive  contribution  to  maintain  the  structural 
information  in  the  resampled  structure,  their  applicability  is  not  completely  assured  in  a  wide  range  of 
grid  sizes.  In  the  last  model,  we  introduce  a  length  ratio  between  the  grid  size  and  the  characteristic  scale 
of  the  original  porous  structure  in  order  to  let  the  model  exert  its  effect  preferably  in  accordance  with  the 
grid  size.  The  model  successfully  represents  the  original  structural  information  in  a  wide  range  of  grid 
sizes  and  improves  the  simulation  accuracy  within  the  given  computational  resources. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  promising  energy  conversion 
systems  because  of  their  high  efficiency  and  fuel  flexibility.  Elec¬ 
trodes  of  the  SOFCs  mainly  consist  of  porous  materials,  and  their 
microstructure  has  a  significant  influence  on  the  electrode  perfor¬ 
mance.  Moreover,  performance  degradation,  which  is  a  significant 
barrier  to  the  commercialization  of  SOFCs,  is  often  attributed  to  the 
morphologic  changes  in  the  porous  structure  [1—7],  such  as 
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sintering  and  redox  of  metal  catalyst  and  carbon  deposition. 
Therefore,  there  is  a  strong  demand  to  clarify  the  effect  of  the 
microstructure  on  the  performance  and  durability  of  the  electrodes. 

Numerical  simulation  is  one  of  the  most  powerful  and  effective 
approaches  to  investigate  the  microstructural  effect  on  the 
performance  [8-14].  Various  phenomena  occurring  inside  the 
electrodes,  such  as  electron/ion  conduction  in  the  solid  phases,  gas 
diffusion  in  the  pore  phase,  and  electrochemical  reaction  at  the 
three-phase  boundary  (TPB),  are  analyzed,  and  the  macroscopic 
performance  of  the  electrodes  is  predicted. 

The  accuracy  of  the  numerical  simulation  basically  depends  on 
two  aspects:  the  electrochemical  reaction  model  and  the  structure 
model.  For  the  modeling  of  the  electrochemical  reaction  around  the 
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TPB,  many  experimental  approaches  have  been  conducted  focusing 
on  the  elementary  electrochemical  reaction  on  the  catalytic  surface 
in  the  electrodes  [15-19].  Some  of  the  findings  were  summarized 
into  numerical  models,  which  are  available  for  the  numerical 
simulation  [20-23].  For  the  structure  model,  on  the  other  hand, 
many  porous  models  have  been  developed  to  mimic  the  complex 
electrode  microstructure  [21,24-26].  For  example,  the  random- 
spherical  packing  model  [24,25]  assumes  that  the  porous  struc¬ 
ture  is  composed  of  two  kinds  of  spheres,  representing  electron- 
conductive  particles  and  ion-conductive  particles,  which  are 
randomly  distributed  in  the  field.  From  the  obtained  structure, 
microstructural  parameters  such  as  phase  connectivity  and  TPB 
density  can  be  obtained.  However,  these  structure  models  require 
assumptions  for  the  structural  configuration,  and  hence,  they 
cannot  always  represent  the  heterogeneous  nature  of  the  complex 
porous  structure. 

Recently,  with  state-of-the-art  observation  techniques  such  as 
scanning  electron  microscopy  combined  with  focused  ion  beam 
milling  (FIB-SEM)  and  X-ray  computed  tomography,  the  three- 
dimensional  (3D)  structure  of  real  porous  electrodes  can  be  ob¬ 
tained  at  nano-scale  resolution  [1,2,4,6,27-31].  A  simulation  tech¬ 
nique  using  the  obtained  3D  structure  has  also  been  recently 
developed  [8,9,12-14].  In  such  3D  analyses,  however,  the  compu¬ 
tational  load  often  becomes  a  critical  issue  when  we  use  the  ob¬ 
tained  voxel-based  structure  for  the  numerical  grid  system.  The 
elemental  image  voxel  in  the  FIB-SEM  observation  is  usually  several 
tens  of  nanometers  in  size;  therefore,  the  total  number  of  grids  tends 
to  reach  to  the  order  of  a  hundred  million,  which  is  unrealistic  for  the 
grid  system  considering  the  available  computational  resources. 
Currently,  it  is  a  common  practice  to  resample  the  original  voxel- 
based  data  set  to  reduce  the  number  of  calculation  grid.  The 
resampling  process,  however,  raises  a  question  as  to  whether  the 
resampled  data  set  is  still  a  good  representation  of  the  original 
microstructure  because  it  inevitably  results  in  the  loss  of  structural 
information  that  is  smaller  than  the  calculation  grid  size  (sub-grid- 
scale  information).  In  our  previous  report  [14],  we  proposed  the  use 
of  a  sub-grid-scale  (SGS)  model  in  order  to  maintain  the  quality  of 
the  original  structural  information  intact  in  the  resampled  structure 
used  for  the  calculation  grid  system.  In  the  proposed  SGS  models,  the 
local  transport  coefficients  in  the  grids  are  evaluated  with  the  local 
information  of  the  porous  structure,  such  as  the  volume  fractions 
and  structural  complexity.  The  effect  of  the  SGS  models  was  inves¬ 
tigated  through  the  diffusion  and  overpotential  analyses.  We  found 
that  if  the  structural  complexity  inside  each  calculation  grid  is 
reasonably  taken  into  account,  the  transport  flux  of  various  species 
in  the  porous  structure  can  be  properly  evaluated;  hence,  a  reason¬ 
able  value  of  the  overpotential  can  be  obtained  even  when  resam¬ 
pling  is  involved  in  the  grid  generation.  However,  the  SGS  models 
examined  in  our  previous  report  did  not  always  exert  the  desirable 
effect  in  a  wide  range  of  grid  sizes:  the  effect  was  not  enough  when 
a  coarser  grid  system  was  applied  or  was  more  than  what  was 
required  in  some  other  cases. 

The  aim  of  this  study  is  to  establish  more  sophisticated  SGS 
models  for  the  3D  numerical  simulation  in  the  porous  structure  and 
expand  its  applicability  limit.  The  local  structural  complexity  inside 
a  grid  is  evaluated  by  the  power  of  the  volume  fraction  of  a  porous 
component  (power  law),  which  has  been  often  used  in  the 
conventional  porous  models  [24,25].  Moreover,  in  order  to  include 
the  effect  of  the  grid  size,  we  introduce  an  adaptive  power  index 
considering  the  relationship  between  the  grid  size  and  the  char¬ 
acteristic  scale  of  the  porous  structure  (power  law  with  an  adaptive 
power  index).  The  effect  of  the  SGS  model  on  the  numerical 
simulation  is  extensively  investigated  through  the  diffusion  anal¬ 
ysis  in  the  porous  structure,  and  their  effect  and  applicability  are 
evaluated. 


2.  FIB-SEM  observation 

2  A.  Sample  preparation 

In  this  study,  a  conventional  Ni-YSZ  cermet  is  used  as  a  sample 
porous  structure  representing  the  SOFC  anode.  The  fabrication 
process  is  as  follows:  first,  NiO  powder  (Wako  Pure  Chemical 
Industries,  Ltd.)  and  YSZ  powder  (Tosoh  Co.)  are  mixed  and  ball- 
milled  with  zirconia  balls  (04.0  mm)  for  24  h  along  with  ethanol 
to  disperse  the  particles.  After  the  milling,  the  ethanol  is  evaporated 
using  a  hot  stirrer,  and  the  resultant  powder  is  pre-sintered  at 
1400  °C  for  5  h.  It  is  subsequently  ground  for  3  h  and  mixed  with 
polyethylene  glycol  to  form  a  slurry.  Finally,  the  anode  slurry  is 
screen-printed  on  a  disk  YSZ  electrolyte  (Tosoh  Co.,  24  mm  in 
diameter,  500  pm  in  thickness)  and  sintered  at  1400  °C  for  5  h.  The 
thickness  of  the  anode  is  around  30  pm.  The  fabricated  NiO-YSZ 
cermet  is  reduced  under  pure  hydrogen  atmosphere  at  1000  °C 
for  1  h  and  then  cooled  down  in  a  reductive  atmosphere.  We 
prepare  anodes  with  three  different  compositions:  Ni:YSZ  =  30:70, 
50:50,  and  70:30  vol.%.  The  reduced  anode  is  infiltrated  with  epoxy 
resin  (Marumoto  Struers  KK)  under  vacuum  condition  so  that  the 
pores  of  the  electrode  can  be  easily  distinguished  during  the  SEM 
observation.  The  infiltrated  sample  is  cut  and  mechanically  pol¬ 
ished  with  sandpaper  and  diamond  paste  to  prepare  them  for  the 
FIB-SEM  observation. 

2.2.  Observation 

The  3D  microstructure  of  the  Ni-YSZ  anodes  is  observed  using 
the  FIB-SEM  system.  An  in-lens  secondary  electron  detector  is  used 
for  the  microstructural  observation  with  an  acceleration  voltage  of 
~2  kV.  Fig.  1(a)  shows  an  example  of  the  obtained  cross-sectional 
images,  in  which  the  white,  gray,  and  black  parts  correspond  to 
Ni,  YSZ,  and  the  pore  phase,  respectively.  In  this  study,  orthogonal 
coordinate  axes  X  and  Yare  embedded  on  the  2D  SEM  image  and  Z 
is  the  proceeding  direction  of  the  FIB  milling. 

We  extract  regions  available  for  the  numerical  analysis  and 
conduct  a  phase  separation  based  on  the  image  brightness.  After 
the  alignment  and  the  phase  separation,  the  sequential  set  of  the 
2D  images  is  lined  up  with  the  actual  increment  in  the  FIB-SEM 
observation  and  the  3D  porous  microstructure  is  reconstructed  in 
a  virtual  field.  Fig.  1(b)  shows  the  reconstructed  microstructure  of 
the  porous  anode  (Ni:YSZ  =  50:50).  We  use  the  commercial  image 
processing  software  Avizo  (Mercury  Computer  Systems,  Inc.)  for 
the  phase  separation,  3D  reconstruction,  and  some  of  the  quanti¬ 
fication  of  the  porous  structure. 

3.  Numerical  model 

We  conduct  a  numerical  simulation  with  the  3D  porous  anodes 
obtained  by  the  FIB-SEM.  In  order  to  evaluate  the  proposed  SGS 
models  in  terms  of  how  precisely  they  can  maintain  the  structural 
information  of  the  original  porous  structure,  an  index  is  required 
to  characterize  both  the  original  porous  structure  and  the 
resampled  structure  used  for  the  calculation  grid  system.  For  this 
purpose,  the  tortuosity  factor  is  the  most  appropriate  micro- 
structural  parameter  because  it  directly  represents  the  complexity 
of  the  porous  structure.  We  apply  different  quantification  methods 
for  the  original  porous  structure  and  the  resampled  structure  as 
follows. 

For  the  original  porous  structure,  we  apply  the  random-walk- 
based  diffusion  simulation  (RW)  [10,32].  This  method  is  special¬ 
ized  for  evaluating  the  tortuosity  factor  from  the  voxel-based 
structure  data  such  as  that  obtained  by  the  FIB-SEM.  Since  the 
requirement  for  the  computational  resources  is  much  smaller  than 
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Fig.  1.  FIB-SEM  observation  of  the  SOFC  porous  anode,  (a)  Example  of  the  cross- 
sectional  images  (Ni:YSZ  =  50:50,  white:  Ni,  gray:  YSZ,  black:  pore),  (b)  Recon¬ 
structed  3D  structure  (Ni:YSZ  =  50:50)  (Green:Ni,  Yellow:YSZ).  (For  interpretation  of 
the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of 
this  article.) 


that  for  the  other  approaches,  resampling  is  not  necessary.  The 
accuracy  of  this  method  has  been  proved  in  our  previous  report 
[29]  by  comparing  the  results  with  those  obtained  by  the  lattice 
Boltzmann  method  (LBM).  Therefore,  the  tortuosity  factor  ob¬ 
tained  by  the  random-walk  simulation  truly  represents  the 
structural  complexity  of  the  original  porous  structure.  In  the  later 
discussion,  we  use  the  values  obtained  from  this  simulation  as 
a  reference. 

For  the  resampled  structure  used  as  the  grid  system,  on  the  other 
hand,  we  apply  the  diffusion  simulation  based  on  the  finite  volume 
method  (FVM)  aided  by  various  SGS  models.  Although  our  in-house 
code  is  capable  of  conducting  the  overpotential  analysis  by  consid¬ 
ering  the  electrochemical  reaction  at  the  TPB  [14],  only  transport 
analysis  is  enough  to  investigate  the  effect  of  the  proposed  SGS 
models.  Fig.  2  shows  the  schematic  diagram  of  the  grid  generation 
for  the  FVM  analysis  using  the  SGS  models.  As  is  often  adopted  in 
most  studies,  the  original  porous  structure  is  resampled  to  have 
a  lower  spatial  resolution.  The  unique  strategy  used  in  our  study  is 
that  the  structural  information  inside  the  grid  is  reflected  onto  the 
resampled  structure  in  the  form  of  effective  transport  coefficients 
such  as  effective  conductivities  and  diffusivities. 

In  the  figure,  the  following  length  ratio  is  introduced  as 
a  parameter  to  generally  evaluate  the  grid  size: 


(i) 


where  L/  is  the  characteristic  scale  of  the  phase  l,  such  as  the 
particle/pore  size,  and  Lgrid  is  the  grid  size  used  in  the  simulation, 
which  is  defined  as  follows: 

Lgrid  =  \/AxAyAz  (2) 

where  Ax,  Ay,  and  Az  are  the  grid  sizes  in  each  direction.  The  value 
of  L//Lgrid  is  interpreted  as  the  number  of  grids  that  are  used  to 
resolve  a  typical  structure  of  a  porous  material,  in  other  words,  one 
particle.  If  the  grid  size  is  smaller  or  L//Lgrid  is  larger,  we  can  expect 
that  the  grid-scale  information  (resampled  structure  without  the 
SGS  model)  represents  the  structural  complexity.  On  the  other 
hand,  if  the  grid  size  is  larger  or  L//Lgrid  is  smaller,  the  proposed  SGS 
models  need  to  exert  their  effects  to  maintain  the  sub-grid-scale 
information  in  the  resampled  structure. 

If  the  tortuosity  factor  obtained  by  the  FVM-based  diffusion 
simulation  is  the  same  as  that  obtained  by  the  random-walk 
simulation,  we  can  judge  that  the  resampled  structure  reasonably 
maintains  the  quality  of  the  original  structure. 


3.1.  Random-walk- based  diffusion  simulation 

The  tortuosity  factor  of  the  original  porous  structure  is  quanti¬ 
fied  by  the  diffusion  simulation  based  on  the  random-walk  process. 
In  the  first  step  of  this  method,  a  large  number  of  imaginary 
particles  are  randomly  distributed  in  the  considered  phase,  i.e.,  the 
Ni,  YSZ,  or  pore  phase.  Then,  each  walker  randomly  chooses  one  of 
the  neighboring  voxels  as  its  possible  location  in  the  next  time  step. 
If  the  selected  neighboring  voxel  is  the  same  as  the  considered 
phase,  the  walker  migrates  to  the  voxel.  If  the  selected  voxel  is 
a  different  phase,  the  walker  stays  at  the  current  voxel  and  waits  for 
the  next  time  step.  In  this  procedure,  neither  absorption  nor 
desorption  on  the  phase  boundaries  is  considered.  While  repeating 
this  process,  the  mean  square  displacement  of  the  random  walkers 
is  calculated  as  follows: 


1  N  r 

<r2(t)>  =  Jj  £  -xn(0)}2+{y„(t)  -yn(  0)}2 

n  =  1 

+  {^n(0  —  Zn(0)}2 


(3) 


where  N  is  the  total  number  of  imaginary  particles,  and  xn(t),  yn(t), 
and  zn(t)  are  the  3D  coordinates  of  the  walker’s  position  at  time  t  for 
the  nth  walker.  ()  indicates  the  average  over  all  the  particles.  Since 
the  mean  square  displacement  r^t)  is  proportional  to  time,  the 
transport  coefficient  in  the  phase  l,  i.e.,  Th  is  related  to  the  time 
derivative  of  r^t)  and  the  volume  fraction  of  the  phase  V/. 


1  d<r2(t)), 

6 


(4) 


The  mean  square  displacement  in  the  porous  media  is  lower 
than  that  obtained  in  a  bulk  material  because  the  movement  of  the 
particles  is  interrupted  at  the  phase  boundaries.  The  degree  of  the 
reduction  is  measured  quantitatively  with  the  tortuosity  factor  % 
which  is  defined  as 


rT  =  v,  d (r2(t))T'k  /d(r2(f))j 
1  rf  vff  dt  /  dt 


(5) 


where  Vff  is  the  volume  fraction  of  the  percolated  cluster. 

The  anisotropic  tortuosity  factor  is  also  evaluated  using  the 
directional  mean  square  displacements  <x2(t)>,  <y2(t)>,  and 
<z2(t)>,  and  a  set  of  similar  equations,  namely,  Eqs.  (3)-(5). 
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Fig.  2.  Schematic  diagram  of  the  grid  generation  for  the  FVM  analysis  using  the  SGS  models. 


3.2.  FVM-based  diffusion  simulation  with  the  SGS  model 
3.2.1.  Governing  equation 

To  evaluate  the  structural  complexity  of  the  resampled  structure 
used  for  the  grid  system,  we  conduct  a  simple  diffusion  simulation 
based  on  the  FVM  using  the  following  equation. 

'oca' V0,)  =  0  (6) 

where  fi  is  an  arbitrary  potential,  which  corresponds  to  the  electric 
potential  or  gas  concentration  in  actual  SOFC  electrodes.  reff,|oca|  is 
the  “local”  transport  coefficient,  which  is  defined  at  every  calcula¬ 
tion  grid  and  includes  the  local  structural  information  evaluated  by 
the  SGS  models.  A  detailed  description  of  the  SGS  models  is 
provided  in  the  next  section.  We  set  the  potential  difference  at  the 
two  facing  edge  surfaces  of  the  calculation  domain  perpendicular  to 
the  Z  axis  and  induce  flux  through  the  calculation  domain  (Fig.  3). 
From  the  amount  of  the  induced  flux,  the  effective  transport 


coefficient  of  the  phase  /,  Tfff,  in  the  whole  calculation  domain  is 
obtained,  and  the  tortuosity  factor  of  the  phase  is  subsequently 
obtained  as  follows: 

7-bulk 

=  Vlfe r  (7) 

1 1 

3.2.2.  Sub-grid-scale  (SGS)  model 

To  consider  the  structural  information  with  a  characteristic  scale 
that  is  smaller  than  the  calculation  grid  size  and  to  evaluate 
j-eff, local ^  we  ^ave  been  deveioping  SGS  models  and  applying  them 
to  the  simulation.  In  this  study,  we  investigate  four  types  of  SGS 
models,  the  first  two  of  which  have  already  been  reported  in  our 
previous  work  [14]. 

3. 2. 2.1.  SGS1:  Volume  conservation.  As  the  simplest  SGS  model,  we 
consider  the  conservation  of  the  phase  volume  in  each  grid  (SGS1 ). 
The  numbers  of  voxels  corresponding  to  the  Ni,  YSZ,  and  pore 
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Fig.  3.  Schematic  picture  of  the  FVM-based  diffusion  simulation. 


phases  are  separately  counted  and  the  volume  fractions  are  ob¬ 
tained.  By  using  the  volume  fractions,  the  local  transport  coeffi¬ 
cients  rf’ local  are  evaluated  in  each  grid  as  follows: 


pe ff, local  _  y local  j^bulk 


(8) 


j-eff, local 

f  (i+  \ii-> 


h(i+2,j,k)  (  1  1  \ 

5  I  reff, local  7-eff, local  I 

all 


(10) 


where  San  is  the  total  surface  area  of  the  grid  interface  and  S/  is  the 
surface  area  of  the  phase  l  on  the  grid  interface.  The  transport 
coefficients  at  the  other  interfaces  are  also  defined  with  a  similar 
formula.  A  more  detailed  description  of  the  SGS2  can  be  found  in 
our  previous  report  [14]. 

3.2.23.  SGS3:  Power  law  with  a  constant  power  index.  Before  3D 
observation  of  the  porous  structure  became  available  with  the  help 
of  FIB-SEM  and  X-ray  CT  technologies,  porous  models  had  been 
used  to  mimic  the  porous  structure.  In  these  models,  the  macro¬ 
scopic  tortuosity  factor  is  often  evaluated  with  the  following 
empirical  formula: 


rf  =  vjrfu[k  (ii) 

where  y  is  the  Bruggeman  factor,  which  represents  the  state  of  the 
tortuous  conduction  paths  inside  the  grids.  A  fully  percolated 
network  with  a  straight  conduction  path  is  represented  as  y=l. 
Since  the  conduction  path  inside  the  porous  structure  is  compli¬ 
cated  and  tortuous,  the  Bruggeman  factor  is  over  unity,  and  usually 
assumed  to  be  1.5-2.0  [25].  In  the  SGS3,  we  use  the  above  macro¬ 
scopic  expression  to  evaluate  the  local  tortuosity  factor  in  each 
calculation  grid  (SGS3): 


where  r|3ulk  is  the  bulk  transport  coefficient  and  V;local  is  the  volume 
fraction  of  the  phase  l  in  the  grid. 

In  the  FVM,  any  transport  between  two  neighboring  grids  occurs 
through  the  interface  they  share.  Since  the  local  transport  coeffi¬ 
cients  of  the  two  grids  are  generally  different,  an  averaged  trans¬ 
port  coefficient  between  the  two  grids  needs  to  be  evaluated  to 
calculate  the  transport  flux  through  the  interface.  Since  an  equally 
spaced  grid  system  is  used  in  this  study,  it  is  evaluated  as 
a  harmonic  average  of  the  effective  transport  coefficients  of  the  two 
neighboring  grids.  The  following  formula  shows  the  averaged 
effective  transport  coefficient  between  the  grid  (i,  j,  k)  and  grid 
(f+1  k). 


j-e  ff, local 

f  (i+ 


( _ ! _ 

1  7-eff, local 

V  IfiJM) 


+ 


1 


peff,  local 

I>0'+1  j,k). 


(9) 


j-eff, local  _  ^yrlocal^ j'bulk  ^2) 

The  volume  fraction,  Vjocal  in  eq.  (12)  is  individually  evaluated  in 
each  calculation  grid,  while  Vi  in  eq.  (11)  is  evaluated  from  the 
whole  structure.  Eq.  (12)  is  equivalent  to  the  following  formula, 
where  the  local  tortuosity  factor  term  is  shown  explicitly: 


reff, local  7-bulk  ^local  A/local V  T  mq) 

*1  =  TlocaT^ 1  ’  Tl  =  [Vl  J  ^) 

Ll 

Here  ijocal  is  the  local  tortuosity  factor  modeled  by  the  SGS3.  The 
local  transport  coefficient  on  a  grid  interface  is  evaluated  with  eq. 
(9).  It  should  be  noted  that  the  SGS3  is  equivalent  to  the  SGS1  when 
the  Bruggeman  factor  is  unity. 


Note  that  the  LHS  of  eq.  (9)  is  defined  at  the  location  of  i+1/2, 
which  corresponds  to  the  interface  of  the  grid  (z,j,  k)  and  grid  (z+1,  j, 
k).  The  transport  coefficients  at  the  other  interfaces  are  also  defined 
with  a  similar  formula. 

3.2.22.  SGS2:  Interfacial  connectivity.  The  coarser  the  grid 
becomes,  the  more  difficult  it  becomes  to  capture  the  complex 
porous  structure  inside  a  grid  only  using  the  volume  fraction  of  the 
porous  components.  The  effect  of  the  structural  complexity  on  the 
transport  phenomena  should  be  considered  in  the  grid.  In  the  SGS 
model  2  (SGS2),  the  surface  fraction  of  the  conductive  volume  is 
evaluated  on  every  grid  interface  and  used  to  represent  the  struc¬ 
tural  complexity.  First,  the  local  transport  coefficient  in  each  grid  is 
evaluated  with  eq.  (8),  and  then  the  value  at  the  interface  is  eval¬ 
uated  using  the  surface  fractions  of  the  conductive  volume  as 
follows. 


3.22.4.  SGS4:  Power  law  with  an  adaptive  power  index.  When  the 
original  porous  structure  is  resampled  to  make  a  grid  system,  the 
heterogeneous  nature  inside  the  grid  is  significantly  different  in 
accordance  with  the  grid  size.  When  the  grid  size  is  larger,  each  grid 
can  contain  several  solid  particles  so  that  relatively  complicated 
transport  path  can  be  created  inside  the  grid.  On  the  other  hand, 
when  the  grid  size  is  smaller,  each  grid  can  contain  only  a  part  of 
the  solid  particles;  hence  the  transport  path  inside  the  grid  is  less 
complicated.  Therefore,  it  may  be  preferable  to  vary  the  effect  of  the 
SGS  models  depending  on  the  grid  size  used  in  the  simulation.  In 
the  SGS4,  we  consider  the  relationship  between  the  grid  size  and 
the  characteristic  scale  of  the  porous  structure.  In  this  sense,  it  is 
reasonable  to  include  the  length  ratio  L//Lgr id  into  the  SGS  model.  In 
this  study,  we  vary  the  Bruggeman  factor  in  the  SGS3  depending  on 
the  ratio  for  further  improvement  of  the  SGS  model. 

The  local  transport  coefficient  and  the  local  tortuosity  factor  are 
defined  as  follows: 
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Table  1  Table  3 

Sample  sizes  and  voxel  sizes.  Average  particle/pore  sizes  measured  by  the  line  intercept  method. 


Sample 

Size  [pm] 

Voxel  size  [nm] 

Sample 

Size  [pm] 

X 

Y 

Z 

X 

Y 

Z 

Ni 

YSZ 

Pore 

Ni:YSZ  =  30:70 

18.8 

11.6 

10.1 

36.2 

36.2 

42.2 

Ni:YSZ  =  30:70 

1.34 

1.74 

1.39 

Ni:YSZ  =  50:50 

18.1 

9.72 

9.97 

34.7 

34.7 

62.3 

Ni:YSZ  =  50:50 

1.87 

1.35 

1.55 

Ni:YSZ  =  70:30 

18.0 

13.2 

10.5 

30.0 

30.0 

52.3 

Ni:YSZ  =  70:30 

2.55 

1.11 

2.19 

where  f  is  a  function  representing  the  dependency  of  y  on  the 
length  ratio.  The  derivation  of  the  function  is  described  in  Section  4. 

4.  Results  and  Discussion 

4  A.  Micro  structural  parameters 

Table  1  shows  the  sample  sizes  and  voxel  sizes  of  the  anode 
samples.  We  successfully  obtained  a  region  larger  than  10  pm3  for 
each  anode  sample,  which  is  large  enough  to  be  representative  of 
the  whole  anode  structure  as  proved  in  our  previous  report  [29]. 
Table  2  shows  the  volume  fraction  of  the  anode  samples.  The  ratios 
of  the  Ni  and  YSZ  volume  fractions,  which  are  Ni:YSZ  =  31.8:68.2, 
53.5:46.5,  and  69.1:30.9,  respectively,  are  close  to  the  intended 
values. 

Table  3  shows  the  average  particle/pore  sizes  measured  by  the 
line  intercept  method  [33,34].  The  values  of  the  solid  phases  have 
a  positive  correlation  with  the  volume  fractions;  a  larger  volume 
fraction  causes  more  sintering  between  the  same  kinds  of  particles. 
The  volume  fraction  of  the  pore  phase  depends  on  that  of  the  Ni 
phase  because  a  large  part  of  the  pore  phase  is  created  by  the 
shrinkage  of  the  NiO  particles  during  the  reduction.  These  particle/ 
pore  diameters  are  used  as  the  characteristic  scale  of  the  porous 
component,  L/. 

Table  4  shows  the  anisotropic  tortuosity  factors  quantified  by 
the  random-walk  simulation.  Since  some  of  the  phases  have  no 
percolating  cluster  in  a  specific  direction,  tortuosity  factors  cannot 
be  evaluated.  Anisotropic  aspects  of  the  samples  are  weaker  than 
those  observed  in  our  previous  report  [10],  and  the  difference  is 
attributed  to  the  larger  sample  sizes,  particularly  in  the  Z  direction. 
In  the  following  discussion,  the  tortuosity  factors  in  the  Z  direction 
are  used  as  the  reference  values. 

4.2.  Diffusion  simulation 

Fig.  4  shows  the  tortuosity  factors  obtained  by  the  FVM-based 
diffusion  simulation  with  the  SGS1,  2,  and  3.  The  tortuosity 
factors  are  normalized  with  the  reference  values  obtained  by  the 
random-walk  simulation.  The  length  ratio  L//Lgr id  is  used  for  the 

Table  2 

Volume  fractions. 


Sample  Volume  fractions  [%] 


Ni 

YSZ 

Pore 

Ni:YSZ  =  30:70 

23.6 

50.6 

25.8 

Ni:YSZ  =  50:50 

34.2 

29.7 

36.1 

Ni:YSZ  =  70:30 

42.3 

18.9 

38.8 

horizontal  axis.  In  every  SGS  model,  the  tortuosity  factor  converges 
to  unity  as  the  grid  size  becomes  smaller.  This  indicates  that  the 
effect  of  the  SGS  models  becomes  smaller  with  finer  grids  and  that 
the  structural  complexity  can  be  represented  simply  by  the  grid- 
scale  information.  In  the  case  of  coarser  grids,  on  the  other  hand, 
only  grid-scale  information  is  not  enough;  hence,  the  sub-grid- 
scale  information  contributes  significantly  to  the  representation 
of  the  structural  complexity.  However,  since  the  method  of  evalu¬ 
ating  sub-grid-scale  information  is  different  among  the  SGS 
models,  the  trends  of  TsGs/'fRW  are  significantly  different. 

First,  in  the  SGS1,  the  values  of  tsgs/'Trw  monotonically  increase 
to  unity  as  the  grid  size  becomes  smaller,  which  is  the  same  trend 
that  is  observed  in  our  previous  report  [14].  The  tortuosity  factors 
are  significantly  underestimated  in  the  case  of  coarser  grids,  indi¬ 
cating  that  only  volume  conservation  is  not  enough  to  represent  the 
structural  complexity  inside  the  grid.  We  need  to  evaluate  a  higher- 
order  quantity  associated  with  the  heterogeneous  nature  of  the 
porous  structure  for  more  sophisticated  SGS  model. 

Secondly,  in  the  SGS2,  the  local  structural  complexity  is  evalu¬ 
ated  with  the  surface  information  of  the  grids.  This  approach 
results  in  drastic  changes  in  the  trends  of  tsgs/^rw-  The  values  of 
^sgs/^rw  approach  unity  even  in  the  case  of  coarser  grids.  Similar 
results  were  also  obtained  in  our  previous  report  [14].  However,  the 
extent  of  the  effect  is  different  from  phase  to  phase.  For  example,  in 
the  YSZ  phase  of  the  Ni:YSZ  =  30:70  anode,  the  tortuosity  factors 
are  reasonably  evaluated  over  a  wide  range  of  grid  sizes.  On  the 
other  hand,  in  the  YSZ  phase  of  the  Ni:YSZ  =  70:30  anode,  the  value 
is  significantly  overestimated  in  the  intermediate  grid  size  region 
(L//Lgrid~5).  These  unsystematic  behaviors  of  the  SGS2  can  be 
attributed  to  the  fact  that  the  relationships  between  the  surface 
information  and  the  inner  information  are  not  straightforward. 
Therefore,  the  SGS2  may  not  always  be  valid  in  a  wide  variety  of 
porous  structures. 

Thirdly,  in  the  SGS3,  the  structural  complexity  inside  a  grid  is 
evaluated  by  the  volume  fraction  and  Bruggeman  factor  y,  whose 
value  is  set  to  1.2, 1.5,  or  1.8  in  this  study.  This  approach  is  to  employ 
the  macroscopic  expression  used  in  the  conventional  porous 
models  (Eq.  (11))  on  the  evaluation  of  the  local  tortuosity  factor  (Eq. 
(12)).  When  the  value  of  y  is  small,  the  structural  complexity  is  still 
underestimated  in  the  coarser  grid  size  region;  the  trends  observed 
in  the  graph  are  similar  to  that  in  the  SGS1.  However,  simply 
increasing  the  value  of  y  is  not  an  effective  solution;  the  structural 
complexity  tends  to  be  overestimated  in  the  case  of  intermediate 
grids  if  a  larger  y  is  applied.  We  found  that  it  was  impossible  to 
properly  represent  the  structural  complexity  over  a  wide  range  of 
grid  sizes  if  a  constant  value  is  applied  for  y. 


Table  4 

Anisotropic  tortuosity  factors  quantified  by  random-walk  simulation. 


Sample 

Ni 

YSZ 

Pore 

X 

Y 

Z 

X 

Y 

Z 

X 

Y 

Z 

Ni:YSZ  =  30:70 

— 

78.30 

— 

1.71 

1.57 

2.57 

6.17 

3.92 

10.82 

Ni:YSZ  =  50:50 

4.91 

3.43 

3.79 

7.93 

5.63 

5.39 

3.48 

3.04 

2.64 

Ni:YSZ  =  70:30 

2.79 

2.41 

1.85 

— 

— 

14.05 

3.26 

2.71 

2.27 

274 


M.  Kishimoto  et  al.  /  Journal  of  Power  Sources  223  (2013)  268— 276 


Fig.  4.  Tortuosity  factors  obtained  by  the  FVM-based  diffusion  simulation  normalized  by  the  reference  tortuosity  factors  obtained  by  the  random-walk  simulation. 


The  SGS1,  2,  and  3  discussed  above  showed  positive  contribu¬ 
tions  to  maintain  the  structural  information  in  the  resampled 
structure  used  as  the  calculation  grid  system.  However,  they  are  not 
always  valid  in  a  wide  range  of  grid  sizes;  underestimation  or 
overestimation  is  often  observed.  This  is  because  the  same  formula 
is  applied  to  all  grid  sizes.  The  heterogeneous  nature  inside  the  grid 
significantly  varies  with  the  grid  sizes  so  that  a  tuning  parameter  is 
required  to  vary  the  effect  of  the  SGS  model  depending  on  the  grid 
size.  This  is  the  motivation  to  consider  the  Bruggeman  factor  y  as  an 
adaptive  power  index  in  the  SGS4. 

In  the  last  model,  we  estimate  the  local  tortuosity  factor  with  an 
adaptive  power  index  y  which  depends  on  the  ratio  of  the  char¬ 
acteristic  length  U  and  the  grid  size  Lgrid-  In  order  to  decide  the 
dependency  of  y  on  the  length  ratio  L//Lgr id  we  regard  y  as  a  fitting 
parameter  and  find  its  appropriate  value  y°  with  which  the 
macroscopic  structural  complexity  is  properly  evaluated;  in  other 
words,  the  value  of  tsgs/trw  becomes  unity.  Fig.  5  shows  the 
required  value  of  y°  for  each  phase  versus  the  length  ratio  L//Lgr id. 
The  calculation  is  conducted  in  all  the  phases  in  all  the  samples 
except  the  Ni  phase  of  the  Ni:YSZ  =  30:70  sample.  Although  there  is 
some  deviation,  the  graph  reveals  an  apparent  trend.  If  the  value  of 
L//Lgrid  is  large  enough  or  the  grid  size  is  small  enough,  y°  converges 
to  unity.  This  indicates  that  the  structural  complexity  of  the  porous 
structure  is  represented  by  the  grid-scale  information  and  that  we 
do  not  need  to  rely  on  the  sub-grid-scale  information  in  conducting 
the  transport  analysis.  On  the  other  hand,  if  the  value  of  L//Lgrid 
becomes  smaller  or  the  grid  size  becomes  larger,  y°  tends  to  steeply 


increase  up  to  around  2.5.  This  result  gives  us  an  important  crite¬ 
rion  in  making  a  calculation  grid  system.  If  the  characteristic  length 
of  the  porous  structure  (e.g.  particle/pore  diameter)  is  represented 
with  more  than  15  grids,  the  structural  information  of  the  porous 
electrodes  can  be  represented  by  the  grid-scale  information. 
Otherwise,  the  SGS  model  should  be  applied  to  the  resampled 
structure  to  maintain  the  quality  of  the  structural  information.  In 
the  anode  overpotential  analysis,  where  we  need  to  consider  the 
transports  in  the  three  phases  at  the  same  time,  the  allowable 
maximum  grid  size  without  the  SGS  model  is  determined  by  the 
structure  that  has  the  smallest  characteristic  scale.  For  the 
Ni:YSZ  =  30:70,  50:50,  and  70:30  anodes  used  in  this  study,  the  grid 


Fig.  5.  Required  Bruggeman  factor  against  the  length  ratio,  L//Lgr id,  and  its  fitted  curve. 
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size  needs  to  be  smaller  than  89,  90,  and  74  nm,  respectively. 
However,  it  is  not  always  possible  to  use  such  small  grid  sizes  when 
we  consider  reasonable  computational  resources,  so  that  the  aid  of 
the  SGS  model  is  significantly  important. 

Considering  the  value  of  the  Bruggeman  factor  must  be  larger 
than  unity,  we  fit  the  obtained  values  of  y°  on  the  following 
exponential  formula. 

y  =  Aex p[  -  A-^- J  +  1  (15) 

\  Lgrid  J 

where  A  is  a  pre-exponential  factor  and  A  is  a  decay  constant.  The 
solid  line  in  Fig.  5  is  the  fitting  curve,  where  A  =  1.44  and  A  =  0.246. 
The  mean  square  deviation  is  0.832.  It  is  worth  noting  the  structure 
data  used  in  the  fitting:  the  volume  fraction  of  the  phase  is  from 
18.9  to  50.6%,  which  covers  structures  both  over  and  below  the 
percolation  threshold  of  ordinary  porous  materials.  The  particle/ 
pore  size  is  from  1.11  to  2.55  pm,  which  is  typical  for  an  SOFC 
electrode  material. 

Fig.  6  shows  the  comparison  of  the  effects  among  the  various 
SGS  models.  All  values  of  tsgs/^rw  obtained  in  all  the  phases  of  all 
the  samples  are  summarized  in  the  graph.  It  can  be  seen  that  the 
SGS4  succeeds  in  making  the  tsgs/^rw  much  closer  to  unity  than  the 
other  models.  Structural  complexity  in  coarser  grid  cases  is  more 
accurately  represented  than  the  SGS1,  and  the  overshoot  seen  in 
the  SGS2  is  moderated.  Although  the  SGS4  does  not  give  perfect 
values  for  in  all  tsgs/^rw  grid  size  cases,  significant  improvement  is 
achieved.  The  attempt  to  include  the  effect  of  L//Lgr id  in  the  SGS 
model  has  a  positive  contribution  to  developing  a  more  sophisti¬ 
cated  SGS  model. 

The  computation  time  required  to  obtain  a  convergent  solution 
with  the  SGS4  is  shown  in  Fig.  7.  All  computation  times  are 
measured  when  using  a  single  core  of  Xeon  W3690  processor.  It 
clearly  shows  the  effectiveness  of  the  SGS  model.  If  a  coarser  grid 
can  be  used  in  the  electrode  simulation  without  losing  accuracy,  the 
electrode  simulation  can  be  applied  to  a  larger-scale  analysis,  or  we 


Fig.  7.  Computation  time  required  to  obtain  a  convergent  solution  when  the  SGS4  is 
applied. 


can  expect  more  accuracy  in  the  simulation  with  limited  compu¬ 
tational  resources. 

5.  Conclusions 

Improving  the  accuracy  of  the  sub-grid-scale  model  used  for  the 
3D  numerical  analysis  in  the  SOFC  electrode  is  accomplished  by 
introducing  an  adaptive  power  index.  Four  types  of  SGS  models  are 
applied  to  the  transport  analysis  to  clarify  their  effect  and  appli¬ 
cability.  Although  the  first  three  models,  namely,  volume  conser¬ 
vation,  interfacial  connectivity,  and  the  power  law  with  a  constant 
power  index  models,  have  a  positive  contribution  to  maintain  the 
structural  information  in  the  resampled  structure  used  for  the  grid 
system,  their  applicability  is  not  completely  assured  in  a  wide  range 
of  grid  sizes.  In  the  last  model,  we  introduce  an  adaptive  power 
index  to  the  power-law-type  model  to  let  the  model  exert  its  effect 
appropriately  in  accordance  with  the  grid  size.  When  the  length 
ratio  between  the  grid  size  and  the  characteristic  scale  of  the 
original  structure  is  adopted  as  an  adaptive  power  index,  the  model 
succeeds  in  representing  the  original  structural  information  in  wide 
a  range  of  grid  sizes  and  reducing  the  computational  cost  without 
losing  simulation  accuracy.  An  accurate  and  effective  SGS  model 
can  promote  a  flexible  application  of  the  numerical  technique  to  the 
analysis  in  porous  materials. 
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Nomenclature 

A  pre-exponential  factor 

L  length  (m) 

r  non-dimensional  displacement 

S  surface  area  (m2) 

V  volume  fraction 

Greek  symbols 
y  Bruggeman  factor 

T  transport  coefficient 

A  decay  constant 

t  tortuosity  factor 

</>  arbitrary  potential 

Subscripts 

i  coordinate  number  in  x  direction 

j  coordinate  number  in  y  direction 

k  coordinate  number  in  z  direction 

l  phase 

Superscripts 

bulk  value  for  bulk  material 

eff  value  for  porous  material 

local  local  value  defined  in  a  calculation  grid 
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